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We study finite-temperature phase structures of liard-core bosons in a two-dimensional optical 
lattice subject to an effective magnetic field by employing the gauged CP^ model. Based on the 
extensive Monte Carlo simulations, we study their phase structures at finite temperatures for several 
values of the magnetic fiux per plaquette of the lattice and mean particle density. Despite the 
presence of the particle number fluctuation, the thermodynamic properties are qualitatively similar 
to those of the frustrated XY model with only the phase as a dynamical variable. This suggests that 
cold atom simulators of the frustrated XY model are available irrespective of the particle filling at 
each site. 
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I. INTRODUCTION 

Ultracold atoms in an optical lattice (OL) have been 
a particularly important field to study a wide range of 
fundamental problems in condensed matter physics 
When an OL is rotated, we can expect versatile cold- 
atom quantum simulators, which can demonstrate vari- 
ous effects caused by a magnetic field such as quantum 
Hall effects This is because the neutral atoms in a 
rotating reference frame experience a Coriolis force of the 
same form as the Lorentz force on charged particles in a 
magnetic field. Recently, two experiments were reported, 
making use of a rotating OL to study quantized vor- 
tices in gaseous Bose- Einstein condensates (BECs) [1, 01 ■ 
Moreover, "synthesis" of a gauge field was realized by 
using a spatially varying Raman coupling between inter- 
nal atomic states to implement the required geometric 
phases opening a door to more wide ranging stud- 

ies associated with an artificial magnetic field that gives 
rise to an orbital motion of neutral atoms. 

The Bose- Hubbard model under an effective magnetic 
field exhibits very interesting physics 0-113, beyond 
the physics of the usual Bose-Hubbard model. These 
properties are inherited from the remarkable structure of 
the energy spectrum for noninteracting problems, i.e., a 
single particle moving on a tight-binding lattice in the 
presence of a uniform magnetic field. Here, the energy 
spectrum depends sensitively on the frustration param- 
eter /, a magnetic flux per plaquette of the lattice, and 
exhibits a fractal structure known as the "Hofstadter but- 
terfly" [IJ]. For a rational value f = p/q (with the inte- 
gers p and q), t here are q bands, and each state is g-fold 
degenerate [25l |26| . A series of theoretical proposals in- 
dicates that it should be possible to implement strong 
gauge fields such as / ~ 1 on an OL p3l-[l3|. In fact, it 
has been experimentally demonstrated recently Q. In 
the strongly interacting regime, where both the parti- 
cle number per site and the magnetic flux per plaquette 
are of order unity, it is theoretically predicted that there 
exist strongly correlated phases representative of the con- 
tinuum quantum Hall states [ll|, [3: 13 • In the weakly- 



interacting condensed phase under the strong magnetic 
fleld, the frustrated Josephson junction arrays provide a 
close analog of that system, and implementations using 
cold atoms have been proposed [13, Ull • 

These studies focused mainly on the ground-state 
properties and the quantum phase transitions at zero 
temperature. However, they are mostly based on the 
exact diagonalization study and are restricted to small 
systems. In this paper, we study the finite-temperature 
phase diagram of hard-core bosons in a two-dimensional 
(2D) OL with an effective gauge field. Our study is rel- 
evant to atoms with low densities and very strong repul- 
sive interactions, e.g., tuned by the Feshbach resonance 
[29l |. In the hard-core limit, the Bose-Hubbard model can 
be mapped to the quantum spin model and described by 
the CP^ (complex projective) operators, which are use- 
ful to construct the path- integral formulation (soj . In 
the high-temperature limit, the quantum Hamiltonian 
reduces to the classical three- component XY model [i.e., 
without nearest-neighbor Sz coupling as shown in Eg. pH)) 
below] frustrated by the gauge field, which is referred to 
as the gauged CP^ model below. This reduction provides 
a practical platform to explore the finite-temperature 
phase diagram of this system and to discuss the detailed 
critical properties of the phase transitions. 

It is known that the 2D frustrated XY model (FXYM), 
which has only two components Sx and Sy, is closely re- 
lated to our model and exhibits very rich properties of 
the phase transition [31I - I4H . For the fully frustrated XY 
model (FFXYM) with / = 1/2, there would be double 
phase transitions, one is associated with the Berezinskii- 
Kosterlitz-Thouless (BKT) transition due to the global 
U(l) symmetry and the other associated with Ising- 
model-like transition due to the Z2 chiral symmetry of 
the ground state. A central issue in the studies of the 
FFXYM has been to clarify how these two distinct types 
of orderings take place [Slj. One possibility is that, even 
at the temperature where the Z2 chirality establishes a 
long-range order at T < Tc, the U(l) phase (XY-spin) 
may remain disordered due to thermally excited, un- 
bound vortices. Then, the orderings of the two vari- 



2 



ables take place at two separate temperatures such that 
Tc > ?BKT- The other possibility is that both order- 
ings of the chirality and the phase take place at the same 
temperature, and the resulting single phase transition is 
neither of the conventional Ising type nor the BKT type 
but follows a new universality class. Many studies have 
discussed this phase transition, but yet do not provide 
conclusive results. The phase transitions for other values 
ofj^, e^ g., / — 1/3, 2/5, etc., were also discussed in Refs. 

Our study is an extension of these studies by incorpo- 
rating the particle number fluctuation (z-componcnt of 
the spin) at each site. It is nontrivial how the additional 
degree of freedom has an influence on the above thermo- 
dynamic properties. To the best of our knowledge, the 
statistical properties of this model have not been con- 
sidered so far. We confine ourselves to the values of the 
frustration parameter f — 0, 1/2, and 2/5 to settle the 
discussion, and make a comparison with the results ob- 
tained by the usual XY model. Although our concern- 
ing model is a classical one (obtained by neglecting the 
quantum fluctuations relevant at low temperatures) , this 
should properly describe the thermodynamic properties 
of the hard-core bosons in an optical lattice at sufficiently 
high temperatures [i^. 

In Secini we describe how to obtain our gauged CP^ 
model, starting from the Bose-Hubbard Hamiltonian 
with the gauge field. The ground-state properties of 
the gauged CP^ model are discussed in Sec. IIIIl In 
Sec. IIVI based on the Monte-Carlo simulations, we study 
the finite-temperature phase structures of the gauged 
CP^ model for the averaged site occupation of hard-core 
bosons being half-filled and non-half-filled. Section Ivl is 
devoted to conclusions and discussion. 



II. MODEL 

A. Hamiltonian for hard-core bosons in an effective 
magnetic field 

We consider a system of iV-bosons put on the sites of 
a 2D square lattice with the size L x L. The 2D Bose- 
Hubbard model subject to a uniform Abelian gauge po- 
tential is described by the Hamiltonian 



H = -- 



t 



U 



2 ^Pi{P^ - 1) 

i 



Here, the operator a^^-* destroys (creates) a boson on the 
lattice site i = (i^jiy). implies a pair of nearest- 

neighbor sites, t, /i, and U{> 0) describe the nearest- 
neighbor tunneling energy, the chemical potential, and 
the onsite repulsion, respectively, pi = a|ai is the num- 
ber operator at i and the Hamiltonian conserves the total 



number of bosons, N = pi. Throughout this work, we 
consider the system without additional trap potentials 
such as a harmonic one. Nevertheless, the results of this 
study can be applied within the local density approxi- 
mation to realistic experimental systems which have an 
additional trapping potential. 

The field Aij describes the imposed gauge potential, 
defined by Aij = ' A ■ dr. All of the physics of the sys- 
tem governed by the Hamiltonian ([Ij is gauge-invariant. 
Hence, its properties depend only on the magnetic fiuxes 
of magnitude B through plaquettes 



$ = 



rfS • B = Aij = BcP = 27r/, (2) 



where a labels the plaquette, and the sum represents the 
directed sum of the gauge fields around that plaquette 
(the discrete version of the line integral), c? is a typical 
lattice spacing and the last equality relates B and /. In 
the following, we use the vector potential in the symmet- 
ric gauge and its form can be written as 



A — (yAx , Ay , Az ) 



B B ^ 

y, — X, 

2 2 



(3) 



This corresponds to a uniform magnetic field B — 
(0,0, -B) in the direction perpendicular to the lattice 
plane. The discrete vector potential is then written as 

'■^ \ T^fix for i = {ix-,iy), j ^{ix-,iy + l)- 

In the following, we take a hard-core limit {U oo), 
where the allowed physical states at i are eigenstates of 
Pi with eigenvalue or 1 and their superpositions. States 
with higher particle number at the same site such as dou- 
ble occupancy are excluded. We introduce a destruction 
operator of the hard-core boson as 0i, which satisfies the 
following mixed canonical-(anti)commutation relations 



0, 



for i ^ j. 



and on the same site. 



= 1. 



(5) 



(6) 



At, 



Thus, the number operator is rewritten as pi 
and its eigenvalue pi is assured to be or 1. Then, the 
Hamiltonian H is written as 



t 



he 



E 



H.c.) 



pY.^\k. (7) 



where p determines the mean density p of hard-core 
bosons per site. 



(8) 
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within the range < p < 1. 

In the hard-core hmit, the Bose-Hubbard model be- 
comes equivalent to a spin- 1/2 quantum magnet, because 
the following relations exist between (pi and the s = 1 /2 
SU(2) spin operator sf'^'^^ 130|, 



1 



st = sf + is" ^61. 



The Hamiltonian ([7]) then becomes 



t 



(9) 



^E^^ (10) 



Here, the conservation of total particle number is in- 
terpreted as the constant magnetization — X^i^i- 
Under Eq.®, the eigenstates have the correspondence 
Iti) = I Pi = 1) a-iid = \pi = 0). This Hamiltonian 
describes a quantum spin-1/2 magnet, experiencing XY 
nearest neighbor spin exchange interactions. These ex- 
change interactions are frustrated due to the gauge field 
Aij. 



B. CP^ variable and path-integral representation 

To study the thermodynamic properties of this sys- 
tem, we evaluate the partition function Z of the grand 
canonical ensemble. 



Z = Tre"'^^'-, /3 



1 



(11) 



We express Z by a path integral which is useful for nu- 
merical calculations. For this purpose, it is convenient 
to introduce a CP^ variable Wi = {'Wii,W2i) G C which 
satisfies the CP^ constraint. 



ki^r + k2zp = 1. 



(12) 



An associated pseudocoherent state \wi) is defined by 

\Wi) =Wu\^i) +W2i\U), (13) 

where \wi) is normalized as {wi\wi) — 1 due to Eq. ([T^. 
and we have generally {wi\w'i) — Wiiw'ii+'w*2'Wi2 = w*w^. 
Let us define the integration measure 



[d^w,]=2 I d'wu / d'w2r5{{w,\w,) - I) (14) 
c Jc 



which satisfies 



j [d^w,]l = 2, y" [d'^w,]w*^Wib = 2 X ij^h = 5ab- (15) 
Then the completeness is expressed as 

/ [d^w,]\w,){w,\ = I t,)(t^ I + I U){U 1 = 1. (16) 



The Hamiltonian can be represented by the CP^ oper- 
ators Wii and W2i which satisfy the bosonic commutation 
relation 

[Wai.Wb.j] = 0, [Wai,wl-] = (5afc%, a,b ^ 1,2, (17) 

and their physical states are restricted as 
X^ w^jWailphys) = Iphys). Here, the correspondence to 
the spin operator is given as 



It*) = wLlvac), |;i)=wTjvac) 



and 



~x,y,z _ 1 / t ',t 



(18) 



(w},,W^Ja^'2''^(wi„U,2*)% 

= w[^W2i, sl^wl^wu, (19) 



where cr^'J''^ are Pauli matrices. Equations (|19p and ^ 
imply the relation (j)i = w\^wii, etc. We also note the 
following relation : 



(20) 



Using these relations and following the standard pro- 
cedure, we can write the partition function Z in the path- 
integral form with the imaginary time r e [0, /3] as 

Z = n/MV(r)]e/o^''^^M, (21) 



where 



4 E «(T)w2i(T)u;2^(T)u;i,(T)e*^- +C.C.) 



+^iY.wUT)wu{T). (22) 



Here, the CP^ operators have been replaced to the com- 
plex numbers by employing the path-integral formula- 
tion. 

To proceed further, we make one simplification by con- 
sidering the finite-T region, such that the r-dependence 
of Wai{T) in the path integral can be i gno red, keeping 
only the zero modes as Wai{T) — ?> Wai [4y|- This cor- 
responds to neglecting the quantum fiuctuations. As a 
result, the problem is reduced to the classical one. Under 
the relations © and (|19l) . one can finally obtain 



(23) 



where Hqpi{w) is the classical version of Eq. ([7]) ex- 
pressed in terms of Wi: 



Hcpiiw) = +C-C.) -A*E^l*^l*' 

= W2iWii. (24) 
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This allows us to study the finite-T phase structure, 
which summarizes the essential properties of the system. 
Besides, the finite-T phase diagram gives a very useful in- 
sight into the phase structure at T = 0; if some ordered 
states are found at finite T, we can naturally expect that 
they persist down to T = 0. Equation (|24|) is referred 
to as the "gauged CP^ model" in the following ^4^. We 
can evaluate the partition function of Eq. (I23p to un- 
derstand the thermodynamic properties of the system, 
using the standard Monte Carlo simulations. From the 
symmetry properties summarized in Appendix [Aj we can 
confine ourselves to0</<l/2 and 1/2 < p < 1 in the 
following argument. 

We note that Eq. ([24|) is viewed as a gauged model of 
three-component normalized classical spin Si [0(3) spin] 
[3 as 

Si = wjawi, Si - = 1, 
[d^w^] = -d^s^ d{s^-s, - 1). (25) 

TT 

C. Relation between the CP^ model and the other 
models 

Our model is an extended version of the FXYM, 

HxY = -JY1 - + )• (26) 

To see this, we rewrite the CP^ variables as 

'~\w2r)- l^sin(V^/2)e^^- j' 

Here, the angle variables have the ranges < ipi < tt, 
< Xii,2i < 27r. The Hamiltonian Eq. ([M)) can be 
written as 

Hqpi = — - ^ sinipi sinipj cos{6i — 6j + Aij) 

-|^cosV^„ (28) 

i 

with 0i = — Xii and the integration measure [d^w^] = 
{in^)-'^ sin ■ipid'ipidXiidX2i. The 0(3) spin of Eq. ^ is 
expressed as 

Si = sin^iCoaOi, = sin V'i sin^^, Si = cos-ipi. (29) 

If we restrict the configuration space with fixed ipi — 7r/2, 
the ground state has (uniform) density p — 1/2, so that 
sf = and all the spins lie in the xy-plane. Then, the 
Hamiltonian reduces to the FXYM [Eq. (HH)]. Here, 
"frustration" refers to the fact that, with / 7^ for any 
plaquette, the angles 9i around this plaquette cannot be 



chosen to maximally satisfy the XY exchange couplings. 
The CP^ model has a site-dependent factor sinipisinipj 
associated with the variation of the particle number at 
each site. Hence, the CP^ model includes the particle 
number fluctuation and goes beyond the XY model based 
on the phase fluctuation only. 

If we take into account the nearest-neighbor repul- 
sive interaction in the Bose-Hubbard model, the hard- 
core constraint yields the gauged spin-half quantum XXZ 
model nil 

-I^Y^st. (30) 

i 

Our model corresponds to the classical version of the 
gauged XXZ model with V = 0, called the XXO (three- 
component XY) model [5^. The XXO model is clearly 
distinct from the XY (two-component XY) model, as the 
spins fluctuate also out of the xy plane. In other words, 
even if the XXO model and the XY model share the same 
form of the Hamiltonian, the associated phase space is 
different. 

In addition, we note that there is a similar model in 
which the spins are random in not only their direction 
but also their magnitudes, known as the "fuzzy" spin 
XY model, ^ 

Hfxy = E ^^^3 cos{0i - 9j). (31) 

While the magnitudes Xi of spins in the fuzzy XY model 
can have any values randomly and continuously site by 
site, the magnitude of 0(3) spin Si at each site in the 
CP^ model is fixed unity, as shown in Eq. (l25t . 

In the following, we compare the thermodynamic prop- 
erties of our gauged CP^ model Eq. and the FXYM 
Eq. (pS)) . The latter has been extensively studied for 
decades [3ll - l4]| . To this end, we have to put the two 
models in the same energy measure by establishing a re- 
lation between J and t. By noting the argument given 
below Eq. (P^)) . let us try to replace the density at site i in 
Eq. (|^ to the average value p, i.e., sin ■0^ sin ijij — sin^ 
with 

sin^ i, = Acof? ^ (^1 - cos^ = 4p(l - p), (32) 

where we have used p ~ (0*0) = (w^wi) = cos^(V'/2). 
This correspondence implies the relation J = tp{l — p). 
We use this relation when the energies in the two mod- 
els (1^51) and (1^51) are compared. Note that this relation 
reflects the particle-hole symmetry and is different from 
the naive replacement J = tp. 
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III. GROUND STATE 

We discuss here the ground-state properties of the CP^ 
model. It is expected that the ground state exhibits sim- 
ilar behaviors to the FXYM, where the magnetic flux 
forms the typical pattern with the q x q unit cell struc- 
tures for / = p/q', the checkerboard pattern emerges for 
/ = 1/2 and the staircase pattern for 1/3 </< 1/2 
[52 - [54| . Since the CP^ model has an additional degree 
of freedom associated with the density fluctuation, it is 
expected that there are some differences from the results 
of the FXYM. 

We calculate numerically the ground state of the 
gauged CP^ model for fixed p by the simulated annealing 
method. Two examples for p — 0.5 and 0.95 are shown 
in Fig. [U where the former corresponds to p = (see 
Appendix |^ while the latter is obtained by adjusting a 
proper value of p for a given /. Figures [1] (a) and[T] (b) 
represent the ground-state energy i^min as a function of 
/. We also plot the energy for the FXYM for compar- 
ison. The shape of the energy curve is non-monotonic 
behavior with respect to /, following the bottom of the 
energy spectrum of the Hofstadter butterfly js^l . 

For p = 0.5 the ground-state energy coincides com- 
pletely with that of the XY model. This is clear because 
ipi is freezed to tt/2 at p — 0.5 and the particle-number 
does not fluctuate. Because sf = there, the ground- 
state properties are not affected at all even for / 7^ 0. As 
p deviates from 0.5, the ground-state energy of the CP^ 
model becomes slightly lower than that of the XY model, 
except for / = 0.5, as shown in Fig. [Ub). 

Figures[Hc) and[T{d) show the distribution of the mean 
density pi and the vorticity mj at the site of the dual 
lattice i, defined by 

Pi = lT.P^ = ^ E - + ^^j)' (33) 

where \9i — 9j + Aij\ < tt. The pattern of rua consti- 
tutes the structure of a q x q unit cell, being similar to 
that of the FXYM. As discussed above, the ground-state 
particle density for p — 0.5 becomes uniform for any val- 
ues of /. For p ^ 0.5, however, the mean density is also 
modulated spatially in accordance with the distribution 
of vorticity except for / = 1/2. This is the reason why 
the ground-state energy of the CP^ model is lower than 
the XY model. 




0.0 0.2 0,4 , 0.6 0.8 1.0 




FIG. 1: The ground-state energy per site of the CP^ model 
Eq. (|28|) as a function of the magnetic flux / for the aver- 
aged density (a) p = 0.5 and (b) 0.95, obtained through the 
simulated annealing. The solid curve denotes the energy of 
the CP^ model without the chemical-potential term, while 
the dashed curve denotes that of the XY model with setting 
J — tp{l — p) for comparison. The two curves overlap com- 
pletely in (a); see the text for the reason, (c) and (d) show the 
distribution over the lattice (L = 12 and 10 for / = 1/2 and 
2/5, respectively) of the mean density (pj) (the upper panels) 
and the vorticity (the lower panels) at each plaquette i 
(the site of the dual lattice) for p = 0.95; (c) f = 1/2 and (d) 
/ = 2/5. 



IV. PHASE STRUCTURES AT FINITE T 

We next turn to the discussion on the finite- 
temperature phase structures of the gauged CP^ model 
and compare the result with the FXYM of Eq. ([26]). As 
described in Sec. H the FFXYM with / = 1/2 may give 
rise to a nontrivial double phase transition |3ll - l4l| as- 
sociated with the Ising transition apart from the BKT 



transition which is normally present at / = 0. In ad- 
dition, there have been some discussions on the phase 
transitions for / = 2/5 and / = 1/3 [43 - |45j : their na- 
ture is dominated by the properties of the domain walls, 
following the Ising-like transition for / = 1/3 and the 
first-order transition for / = 2/5. Our primary interest 
is to see how an additional degree of freedom, associated 
with the particle number fiuctuation at each site, mod- 
ifies these properties. We made multicanonical Monte 
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Carlo simulations [55| to calculate statistical averages of 
some quantities described below. 

We study the thermodynamic properties by calculating 
the specific heat defined by 



1 



(34) 




This value can be useful to locate the first- and second- 
order phase transition. In addition, to study the BKT 
transition, we calculate the in-plane susceptibility defined 
as X = d{(f>)/d{(3h)\h~^Q with the site-average of the hard- 
core boson field (j) — ipi/V , explicitly written as 



(35) 



Here, an auxiliary term ~h^^ (pi was introduced in Eq. 
([M)) to derive Eq. ([35]). The second term of the right- 
hand side of Eq. ([55)1 vanishes because of the global 
U(f ) symmetry. This value grows as with the 

system size L when the correlation function obeys the 
power-law behavior {(t>*(t)j) (with rj < 2) due to 

the presence of the quasi-long range order. On the other 
hand, it remains finite for L — ^ oo when the correlation 
decays exponentially as oc e"™''. The value 77 

is an exponent for the in-plane correlations below the 
BKT transition temperature Tbkt, being dependent on 
the temperature. We can calculate the exponent, as a 
function of temperature, by fitting the susceptibility for 
several system sizes L to the above expression for each 
temperature. For the conventional XY model the critical 
temperature Tbkt can be estimated when x grows as 
L^/"*, i.e., 7? = 1/4 

Furthermore, we study the helicity modulus which is 
directly connected to the superfluid density. The helicity 
modulus T is a measure of the resistance to an infinites- 
imal spin twist across the system along one coordi- 
nate. More precisely, it is defined through the change 
of the total free energy F with respect to an infinitesi- 
mal twist on the spin configuration along, say, the x-axis 
Qi — Oi' Oi — Oil + 66, where i' is the nearest-neighbor 
site of the i-site along the a;-direction and 69 — Ad/L. 
One readily finds 



AF = T{6ef + 0{{69)'^), 



and T can be given as 



1 



(36) 



(37) 



where is the a;-bond part of the Hamiltonian at AO = 
and Ix is the total current in the x-direction. For the 
CP^ model, they are written as 

ffi = ^ sini/'i sm-ipj cos{Ot - Oj + Aij), 
Ix ~ —-^ ^ sintpi sin^j sin(0i — Oj + Aij). (38) 




— /=0 

— /= 1/2 
— /= 2/5 
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fit 



FIG. 2: The mean density p and the particle number fluc- 
tuation Ap as a function of at the high-temperature or 
zero-hopping limit /3t — > 0; the behavior is thus independent 
of/. 



According to the renormalization-group theory, the helic- 
ity modulus for the conventional XY model in an infinite 
system jumps from zero to the finite value (2/7r)fcB2BKT 
at the critical temperature T — Tbkt [56]. Therefore, a 
rough estimate of the critical temperature at a finite sys- 
tem could be obtained simply by locating the intersection 
of T as a function of T and the straight line T = 2k-BT/'K. 

In the FFXYM, the jump size has been suggested to 
be lager than the universal value (3ll - [33 , l4l| . As a more 
accurate method to estimate both Tbkt and the jump 
size in the helicity modulus, an useful finite-size-scaling 
expression at T = Tbkt is known as [57] 



T(i,r = Ti 



BKT) 



-fcRT, 



BKT 



1 



1 



1 



2 In L 



(39) 



with Tgj^rp and c being fitting parameters. Tgj^rp is re- 
lated to the jump size, being equal to Tbkt in the case 
of the usual XY model. By making a fit of the numer- 
ical data to Eq. (j39p at various temperatures, one can 
estimate the transition temperature Tbkt as well as the 
critical exponent from the jump size as r] — TBKT/4Tgp,rp 

Ellsi. 



A. Density fluctuation 

Before discussing the detailed thermodynamic proper- 
ties, let us see the magnitude of the spatial fiuctuation 
of the particle density, which is the important difference 
between the gauged CP^ model and the FXYM. The par- 
ticle density fiuctuation is defined as 



(40) 



which is zero for the XY model. In the high-temperature 
limit /3t — > 0, one can calculate exactly the partition 
function Z = [(e'^'' - 1)//3m]^, and thus 



/3m 



_ l + (-l + /3/i)& 
^~ /?/i(-l + e/3A') ' 

-2 + (2 - /3V^)e^^ 
^VFTT^^^^) ' 



(41) 
(42) 
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FIG. 3: The particle number fluctuation Ap as a function 
of pt for (a) /i = and (b) /i = 0.7. The solid, dashed, and 
dotted curves correspond to / = 0, 1/2, and 2/5, respectively. 



B. The case of half-filling (^t = 0) 

Here, we consider the case of half-filling p = 0.5 by set- 
ting p — 0. Then, the density distribution is completely 
uniform in the ground state, as seen in the complete over- 
lap of the ground-state energy in Fig. [Ija), and thus the 
CP'^ model reproduces the ground state of the XY model. 
Hovifever, one has to take into account the density fluc- 
tuation at finite temperatures. This additional degree 
of freedom makes the finite-temperature phase diagram 
and the nature of the phase transition nontrivial. Here 
we focus on the case / = 0, 1/2 and 2/5, each of which 
has been known to give rise to quite different phase tran- 
sitions in the FXYM. 



1. f=0 

First, we consider the situation of a zero magnetic field 
f — 0, which is useful to confirm the accuracy of our nu- 
merical computation. Then, the model is equivalent to 
the XXO (three-component XY) model studied in Ref. 
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which are shown in Fig. [21 As one can see below [Fig. 
[9Ka)], the mean density p is weakly dependent on /St for 
/3p 7^ 0. The density fluctuation Ap has a maximum 
at /3p — 0, corresponding to half-filling p = 1/2, and 
decreases for (3fj, ±oo. Since this Ap gives an upper 
limit of the expected particle number fluctuation, one 
can see that influence of the particle number fluctuation 
is less than 20% in a temperature range of our interest. 

As seen in the ground state (Sec. IIH[) . the density 
becomes uniform for — 0. Thus, Ap should go to zero as 
(3t ^ ooioT p, = 0. On the other hands, for ^ ^ 0, it must 
remains finite because the ground state possesses spatial 
density modulation. Figure [3] shows the /3t-dependance , ^ 
of Ap for /J, = and 0.7 for several values of /. For v^J 0.6 
p — 0, Ap approaches to zero a.s pt for any values 
of /. For p = 0.7, on the other hand, Ap approaches to 
zero a.s pt only for / = 1/2, but remains finite for 
the other values of /. This behavior is consistent with 
the ground-state property shown in Fig. [TJ 
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FIG. 4: Several thermodynamic quantities for /=0: (a) The 
specific heat C, (b) the in-plane susceptibility Xi and (c) the 
helicity modulus T for the CP^ model with p = and / = 
as a function of jSt. The system size for each curve is L = 12, 
24, 36, 48. An additional curve T = S/Trpt in (c) indicates 
the universal jump of a BKT transition. 



[50| ■ The phase transition of this model has been found 
to be consistent with the BKT theory. The specific heat 
C has very small finite-size effects. The in-plane suscep- 
tibility X is a strongly increasing function of the system 
size L for low temperatures, while all the data fall on the 
same curve for high temperatures. These facts indicate 
the absence of second-order transition and the possibility 
of the BKT type transition involving the power-law decay 
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of the correlation function. The transition temperature 
of the XXO model was obtained as Tbkt — 0.699J/fcB 
[soj . which is lower than the usual (two-component) XY 
model Tbkt = 0.898 J/fce [Hg"]. This is naturally under- 
stood due to the difference of the degree of freedom of 
these models. 

Our numerical result is shown in Fig. |4l The size de- 
pendance of C is actually small. We estimate the BKT 
transition temperature Tbkt by two methods, one by 
using the data of x and the other by using T. We plot 
x/L'^^'^ for some values of L as a function of I3t by as- 
suming rj — 1 /4, in which the crossing point of the curves 
gives Tbkt- From T, we take the temperatures that 
correspond to the crossing points of T(T) and the line 
T = (2/7r)fcBTBKT related to the universal jump value 
for various values of L, interpolating them to T — > cxd. 
Both these methods give the same T « 0.702 J/Zcb, which 
is consistent with the result of Ref. [50^ and confirms the 
accuracy of our numerical computations. 



2. f=l/2 

We next consider the case with full frustration / = 1/2. 
If there is a continuous phase transition, various quanti- 
ties should exhibit a singular behavior near the transition 
temperature T^- Figure [Hla) represents the specific heat 
C as a function of where C exhibits a peak struc- 
ture as a function of fit and the height of the peak in- 
creases with increasing L. This suggests the occurrence 
of a second-order phase transition. Concurrently, x and 
T grow from zero with increasing /3t, which is a signa- 
ture of the emergence of superfluid order. Hence, the 
qualitative feature of the phase transition is similar to 
the FFXYM. 

It is important to clarify whether the nature of the 
phase transition is consistent with those observed in the 
analysis of the FFXYM. In the FFXYM, two separate 
phase transitions may occur, corresponding to the break- 
ing of the Z2 chirality and the U(l) symmetry (3ll - l4]| . It 
is still inconclusive whether the former obeys the univer- 
sality class of the Ising transition and the latter is subject 
to the BKT mechanism with non-universal jump of the 
helicity modulus. 

First, we focus on the data of the specific heat to clarify 
the Z2-related phase transition. To determine the critical 
temperature and the critical exponents, we use finite- 
size scaling analysis [g^], where the specific heat C is 
expressed as 



(a) 



C{L,t) = L'^'^ct^iL'^H), 



\lvl 



(43) 



where i = (T — Tc)/Tc is a reduced temperature, a the 
standard exponent for the specific heat, and v the expo- 
nent for the divergence of the correlation length. Since 
the scaling functions (j> should depend on a single vari- 
able, we can make all the data for each system size L fall 
on the same curve by appropriately adjusting the values 
of the critical exponents a, v and Tc. For a finite lattice 
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FIG. 5: Several thermodynamic quantities for f—1/2: (a) 
The specific heat C, (b) the in-plane susceptibility Xj and (c) 
the helicity modulus T for the CP^ model with ^ = and 
/ = 1/2 as a function of The system size for each curve 
is L = 12, 24, 36, 48, 60. An additional curve T = S/nfit in 
(c) indicates the universal jump of a BKT transition. 



the peak in the specific heat scales with system size as 
Cmax oc T"/'^ and occurs at the temperature where the 
scaling function (f){L^^'^i) is maximum, defining the finite- 
lattice transition temperature Tc(L) — Tc+constxL^^^'^ . 

The obtained scaling functions are plotted in Fig. [5] 
for both the FFXYM and the gauged CP^ model. For 
the FFXYM, we obtain T^ = Q.454J/fcB, = 0.873 and 
a = 0.383 from Fig. ^a.), which is consistent with the 
hyperscaling relation dv = 2 — a {d is the dimension 
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FIG. 6: The left panels show the power-law scaling collapse 
of the specific heat data for (a) FFXYM and (c) gauged CP^ 
model. The right panels show the power-law and logarithmic 
fitting of the specific heat peak with respect to sizes L for 
(b) FFXYM and (d) the gauged CP^ model. The power-law 
fitted value is a/v — 0.439 for (a) and a/u — 0.397 for (c). 



number) . Our scaling analyses support the non-Ising ex- 
ponent v < 1, which is consistent with some of the lit- 
erature [H-lij, HI]. However, it should be noted that 
there have been several claims that this non-Ising ex- 
ponent is caused by the artifact of the finite-size effect 
and the exponent in the infinite system may be the Ising 
one 1/ = 1 [H, m, Hlj. Since our calculation does not 
have enough system size to resolve this problem and we 
cannot distinguish the data of Fig. iGja) as a power-law 
fitting or logarithmic fitting [see Fig. IHb)], we shall not 
go to discuss the details of this issue. Our main claim in 
this work is that the similar behavior also occurs for the 
gauged CP^ model. For this model, we also extract the 
critical exponents from the same analysis for Fig. |5] (c) 
and (d) as = 0.369J/fcB, v = 0.878 and a = 0.349. 
The transition temperature becomes lower than that for 
the FFXYM, while the critical exponents are similar. 

Next, we study the U(l)-related phase transition by 
employing the same analysis for the / = case. Some 
studies revealed that the jump size of T at T = Tbkt may 
be non- universal for the FFXYM [siMsllilll . Here, we do 
not assume rj = 1/4 and evaluate Tbkt with two different 
methods: (i) Using the x^-fit of the scaling relation Eq. 
((39)l . we evaluate Tbkt and the jump size Tbkt ^bkt, 
which gives the exponent 77 — Tbkt /4Tgj^rp . (ii) We 
plot x/L'^~^ as a function of /3t for several system sizes 
L and take the temperature at the crossing point. We 
make this analysis by varying 77 around 1/4 and search 
the value of 77 that gives the same Tbkt obtained in the 
analysis (i). The summary of this analysis is shown in 
Fig. [71 For the FFXYM, the scaling analysis (i) alone 
gives Tbkt — 0.437J/fcB and 77 = 0.2. This is consistent 
with the previous literature, where Tbkt — 0.437 J/fcB is 
slightly lower than and the BKT jump is non-universal 
[3ll - l34l . l4l| . However, the crossing point obtained by the 
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FIG. 7: Data of the BKT transition of the FFXYM [(a) and 
(c)] and the gauged CP^ model for / = 1/2 [(b) and (d)]. In 
(a) and (b), we show Tbkt obtained by method (i) (see the 
text) by the horizontal line; (a) T = 0.437 J/fce and (b) T = 
0.363J/fcB . In addition, the temperature corresponding to the 
crossing point in the method (ii) is plotted as a function of 77. 
Figures (c) and (d) show the best fit of the scaling relation 
Eq. ([39} at the corresponding Tbkt, which is obtained by 
finding the minimum of x^-fit error shown in the inset. 



analysis (ii) is preferable to 77 « 0.25, as shown in Fig. 
[7fa), which suggests the same universality of the con- 
ventional BKT transition. This usual BKT behavior for 
the FFXYM was also suggested by Olsson [35*]. Simi- 
lar behavior is also found for the gauged CP^ model as 
Tbkt = 0.363 J//cb and 77 « 0.22 for the analysis (i) and 
ri ~ 0.26 for the analysis (ii). 



3. f=2/5 

The thermodynamic properties for / = 2/5 would ap- 
pear to be similar to the / = 1/2 situation as seen in Fig. 
[51 However, the nature of the phase transition is very 
different. In the FXYM, several works indicated that 
the transition is associated with first-order type [i^, |4^ . 
Li and Teitel observed hysteresis of the internal energy 
when the temperature was cycled around the transition 
and used this as an argument for a first-order transition 
[4^ . Denniston and Tang pointed out that the compli- 
cated branching structure of domain walls is similar to 
the q > 5 Pott's models where the first-order transition 
occurs 45]. The most direct indication of a first-order 
transition is the presence of a free energy barrier be- 
tween the ordered and disordered states which diverges 
as the system size increases [i^ . Since there is no diverg- 
ing characteristic length to which the linear dimension L 
could be compared at a first-order transition, one finds 
that it is simply the volume that controls the size 
effects. 

In the gauged CP^ model, Fig. ^a) clearly shows the 
rapid growth of the peak of C. The inset of Fig. ^h) 
shows the peak values of C as a function of L^. The linear 
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FIG. 8: Several thermodynamic quantities for the CP^ model 
with fi — and / = 2/5 as a function of pt. (a) The specific 
heat C, (b) the in-plane susceptibility, and (c) the helicity 
modulus T. The system size for each curve is L = 20, 30, 40, 
50. (d) represents specific heat vs 



fit clearly shows the expected first-order scaling behavior. 
From the positions of the peaks as a function of L, we 
obtain Tc = 0.193 J/Zcb, which is again slightly lower than 
that of the FXYM = 0.2127 J/Zcb [H- The growing x 
and T at low temperature certainly provides the emer- 
gence of the superfluid order. Due to the presence of the 
first-order transition, it is difficult to explicitly discuss 
the properties of the BKT transition. 



C. The case of non half- filling (/i 0) 

Finally, we show the similar data with the previous 
subsection but for fi ^ 0, where the particle occupation 
at each site is not half-filling. In terms of the pseudospin 
Hamiltonian pop . this situation corresponds to applying 
a longitudinal magnetic field, as seen in the second term 
ofEq. (fTUl) . 

We make the similar analysis as in the previous cal- 
culation for /i = 0.7, where the averaged density is 
about p « 0.7-0.8, which depends weakly on / and f3t. 
The thermodynamic quantities behave similarly to those 
found in the half-filling case; an example for fi — 0.7 and 
/ = 1/2 is shown in Fig. |9l We find no qualitative dif- 
ference in the phase transition for each / from the /i = 
case. We show in Table [T] the obtained critical exponents 
u, a and ?7(rBKT) and the critical temperatures Tc and 
3bkt associated with the Z2 and U(l) symmetry break- 
ing, respectively. One can see that the transition tem- 
perature is further reduced from the fi = case. Also, 
the critical exponents are modified from the values of 
/i = 0. This is naturally understood as follows. As one 
increases the average magnitude of the XY spin com- 
ponent [{sf y + (sf )^]^/^ decreases because \sf \ increases 
(For example, the limit fi ±00 implies Sz = ±1). If 
the length of XY spin is fixed, the model should be in 
the same universality class as the FXYM. However, our 
study for /X = and / = 1/2 exhibits that the fluctua- 
tions of s^'^ give rise to critical exponents different from 
those of the FFXYM. Therefore fluctuations around the 
XY spins of different length may certainly produce dif- 
ferent critical exponents. 



V. CONCLUSION AND DISCUSSION 

We study the finite-temperature phase structures of 
hard-core bosons in a two-dimensional optical lattice 
subject to an effective magnetic field by employing the 
gauged CP^ model. Based on the multicanonical Monte 
Carlo simulations, we study their phase structures at fi- 
nite temperatures for several values of the magnetic flux 
per plaquette of the lattice and mean particle density. A 
summary of this work is listed in Table ID Also, the mag- 
nitudes of the particle density fluctuations are measured 
to be less than 20%. They cause the shift of transition 
temperatures Tc and Tbkt, which are slightly decreased 
from those of the FXYM, and also the shift of critical 
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TABLE I: List of the transition temperatures and some critical exponents obtained in this work. The transition temperature 
is measured by using J = tp{l — p); for /i 7^ 0, p is used at the corresponding transition temperature. In 77(rBKT), we represent 
two values obtained by the method (i) and (ii) in Sec. IIV B 21 For comparison, the corresponding values for the 2D Ising model 
are v — 1 and a = 0, which implies the logarithmic divergence of the specific heat. 

Tc a Tbkt ??(rBKT) 

CP' model, / = 0, ^ = = ~ = 0.702(1) J/fce 025 

CP^ model, / = 1/2, ^ = 0.369(4) J/ks 0.878(5) 0.349(9) 0.363(3) J/fce (i) 0.22, (n) 0.26 

CP^ model, / = 2/5, ^ = 0.193(1) J/ks _ _ _ _ 

CP^ model, / = 0, ^ = 0.7 — — — 0.672(1) J/fce 0.25 

CP^ model, / = 1/2, /i = 0.7 0.347(1) J/ks 0.781(25) 0.381(7) 0.334(4) J/fca (i) 0.19, (ii) 0.24 

CP^ model, / = 2/5, n = 0.7 0.168(1) J/ks _ _ _ „ 

2D XY model — — — 0.898(1) J/fca 0.25 

2D FFXYM 0.454(1) J/ks 0.873(3) 0.383(34) 0.437(3) J/ks (i) 0.2, (ii) 0.25 



exponents a and v. However these fluctuations do not 
modify the global phase structure and the critical prop- 
erties of FXYM. 

The regime described by the XY model (Josephson 
junction array) can be realized when the mean parti- 
cle number at each site is very large pi 3> 1 [M, HII, 
where the particle number fluctuation becomes negligi- 

— 1/2 

ble as ~ Pj . The BKT transition in such a regime 
was observed by Schweikhard et al. [63 |. The impor- 
tant message of this work is that, even though the strong 
particle number fluctuation becomes remarkable due to 
the small site occupation pi 1, one can expect similar 
thermal phase transitions seen in the FXYM. The recent 
experimental demonstration on generating an effective 
magnetic field in an optical lattice Q opens the door to 
explore the rich finite-temperature phase diagram of this 
system. 

Acknowledgments 

The authors thank Shu Tanaka for useful discus- 
sions. One of the authors (K.K.) is supported in part 
by a Grant-in-Aid for Scientific Research (Grant No. 
21740267) from MEXT, Japan. 

Appendix A: Symmetry of the gauged CP^ model 

We summarize the symmetry properties of the gauged 
CP^ model Eq. (|24)) . from which one can get some useful 
information to understand the results. 

The model equation ([M)) or has the following sym- 
metry properties: 

1. The model has global U(l) symmetry; it is invariant 
under the change 4>i cf)'^ — e^^(f)i. 

2. The model also has local gauge symmetry. If we 
change the gauge A A + Vx, then the Hamil- 
tonian remains unchanged if the boson picks up 
a phase change as Aij Aij + {xj — Xi) ^^'^ 

— ^ e*-^*0i. In terms of the pseudospin Si for the 



i-th site, this corresponds to a rotation with the 
angle Xi in the xy plane, sf — )> e^*-^'s^. Because 
the choice Xi = T^ixiy gives rise to a shift of fluxes 
/ ,/ + 1, the system is periodic in / with the 
period 1. 

3. The Hamiltonian Eq. (l24l) is invariant under the 
change = corresponding to the time re- 
versal operation 0^ 0' = 0* (A^^ A^^ = -\ai)- 
Since the partition function is not affected by this 
transformation, one can show that the internal en- 
ergy E = (H) = —din Z/d/3 and the mean number 
density p = {N)/L^ = L-^dlnZ/d{Pp) have the 
following properties 

E{(i,t,iiJ)=E{p,t,p,~f), (Al) 
p{p,t,pj) = p{(3,t,p,-f). (A2) 

The symmetric form of the ground state energy 
shown in Fig. [T]can be understood as follows. The 
plotted energy Erain in Fig. [T]is the first term of Eq. 
((M)) E = E^in-ppL'^- From Eq. (K2\i . one can see 
pip{f),(3,tj) = p{pi-f),P,t,-f). For the fixed 
p, Eq. jJT]) gives 

En^init, p(/), /) = i?min(i, -/) 

= i?,ni„(i,M(l -/),!-/). (A3) 

The last equality is due to the invariance for / — > 
/ + !• 

4. Let us consider the transformation ij^i ^ ip[ = -k — 
ipi. Then, Eq. ^ yields i/cpi(Ai) ^ H'^p,{p) = 
Hcp^{~p) and p = (cos^(i/'/2)) -J> p'{p) = 1 - 
p{—p) (or equivalently sf — s'^^ — — sf). This re- 
flects the particle-hole symmetry, where p' should 
be interpreted as the hole density if p represents the 
particle density. Therefore, especially for /i = 0, we 
obtain the relation 

p(/3, t, 0, /) = p'ip, 0, /) = 1 - p(/3, 0, /) - i. (A4) 
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FIG. 9: Several thermodynamic quantities for the CP^ model 
with fi = 0.7 and / = 1/2 as a function of pt. (a) mean 
particle number p (b) The specific heat C, (c) the in-plane 
susceptibility Xi and (d) the helicity modulus T. The system 
size for each curve is L = 12, 24, 36, 48. 
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